{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"import numpy as np\n",
"\n",
"from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n",
"plt.rcParams[\"figure.figsize\"] = (6,5) # set default size for all figures\n",
"\n",
"from scipy.integrate import quad\n",
"from scipy.optimize import brentq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Mean field and RPA response for the one dimensional Hubbard model.\n",
"\n",
"Comparing numerical calculation with rank 4 tensor and analytical results from Fazekas (Ch. 7.4)\n",
"\n",
"Author: Hugo U. R. Strand (2018)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the one dimensional Hubbard model with nearest neighbour hopping\n",
"\n",
"$$\n",
"H = -t \\sum_{i\\sigma} (c^\\dagger_{i\\sigma} c_{i+1\\sigma} + c^\\dagger_{i+1\\sigma} c_{i\\sigma} ) + U \\sum_{i} \\hat{n}_{i\\uparrow} \\hat{n}_{i\\downarrow}\n",
"$$\n",
"\n",
"The non-interacting momentum space dispersion is $\\epsilon(k) = -2t \\cos k$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mean-field decoupling\n",
"\n",
"$$\n",
"\\hat{n}_\\sigma = \\delta\\hat{n}_\\sigma + n_\\sigma\n",
"$$\n",
"\n",
"$$\n",
"n_\\sigma = \\langle \\hat{n}_\\sigma \\rangle = \\frac{n}{2} \\pm m\n",
"$$\n",
"\n",
"$$\n",
"\\hat{n}_{\\uparrow} \\hat{n}_{\\downarrow}\n",
"=\n",
"\\delta\\hat{n}_\\uparrow\n",
"\\delta\\hat{n}_\\downarrow\n",
"+\n",
"\\big(\\frac{n}{2} - m\\big) \\hat{n}_\\uparrow\n",
"+\n",
"\\big(\\frac{n}{2} + m\\big) \\hat{n}_\\downarrow\n",
"- \\frac{n^2}{4} + m^2\n",
"$$\n",
"\n",
"neglecting fluctuations\n",
"\n",
"$$\n",
"H_{MF} =\n",
"\\sum_{k\\sigma} \\Big[ \\epsilon(k) + \\frac{U}{2}n \\mp Um \\Big] c^\\dagger_{k\\sigma} c_{k\\sigma}\n",
"- U \\Big( \\frac{n^2}{4} - m^2 \\Big)\n",
"$$\n",
"\n",
"\n",
"$$\n",
"E[n, m] =\n",
"\\sum_\\sigma\n",
"\\int_{-\\pi}^{\\pi} \\frac{dk}{2\\pi}\n",
"\\tilde{\\epsilon}_\\sigma(k) f(\\tilde{\\epsilon}_\\sigma(k))\n",
"- U \\Big( \\frac{n^2}{4} - m^2 \\Big)\n",
"$$\n",
"\n",
"with $\\tilde{\\epsilon}_\\sigma(k) = \\epsilon(k) + \\frac{U}{2}n \\mp Um$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def get_total_energy_mf_ref(t, beta, U, mu, n, m):\n",
"\n",
" def disp(k, U, n, m, s):\n",
" return -2*t * np.cos(k) - s*U*m + 0.5*U*n\n",
" \n",
" def integrand(k, U, mu, n, m, s):\n",
" e = disp(k, U, n, m, s)\n",
" f = 1./( np.exp(beta * (e - mu)) + 1 )\n",
" return e * f\n",
"\n",
" E_kin_up, err = quad(integrand, -np.pi, np.pi, args=(U, mu, n, m, +1.))\n",
" E_kin_do, err = quad(integrand, -np.pi, np.pi, args=(U, mu, n, m, -1.))\n",
" E_kin = E_kin_up + E_kin_do\n",
" k_vol = 2*np.pi \n",
" E_kin *= 1. / k_vol\n",
"\n",
" E_tot = E_kin - U*(0.25*n**2 - m**2)\n",
" \n",
" return E_tot"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = 1.0\n",
"beta = 5.0\n",
"n = 1.0\n",
"\n",
"U_vec = np.array([0., 1., 2., 3., 4., 5., 6., 7., 8.])\n",
"m_vec = np.linspace(-0.5, 0.5, num=100)\n",
"\n",
"for U in U_vec:\n",
" mu = 0.5*U \n",
" E_vec = np.array(\n",
" [ get_total_energy_mf_ref(t, beta, U, mu, n, m) \n",
" for m in m_vec ])\n",
" plt.plot(m_vec, E_vec, '-', lw=0.5, label='U=%2.2f' % U)\n",
"\n",
"plt.title(r'$\\beta=%2.2f$, $n=%2.2f$, $t=%2.2f$' % (beta, n, t))\n",
"plt.xlabel(r'$m$')\n",
"plt.ylabel(r'$E[m]$')\n",
"plt.legend(loc='upper right')\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Analytic susceptibility\n",
"\n",
"Density of states\n",
"\n",
"$$\n",
"\\rho(\\epsilon) = \\frac{1}{2t \\pi} \\frac{1}{ \\sqrt{1 - \\epsilon^2/(2t)^2 }}\n",
"$$\n",
"\n",
"$f(\\epsilon) = [e^{\\beta \\epsilon} + 1]^{-1}$, \n",
"\n",
"$$\\partial_\\epsilon f(\\epsilon) = -\\frac{\\beta}{4 \\cosh^2(\\beta \\epsilon/2) }$$\n",
"\n",
"$$\n",
"\\chi_0(k=0) =\n",
"- \\int_{-2t}^{2t} d\\epsilon \\rho(\\epsilon) \\partial_\\epsilon f(\\epsilon) \n",
"$$\n",
"\n",
"Instability rule\n",
"\n",
"$$\n",
"1 - U \\chi_0(k=0) = 0\n",
"$$\n",
"\n",
"$$\n",
"\\chi = \\frac{1}{2} \\frac{\\chi_0}{1 - U \\chi_0}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def get_density_of_states(eps, t):\n",
" return 1./np.sqrt(1. - (eps/(2*t))**2) / (2*t * np.pi)\n",
"\n",
"def fermi_distribution(beta, eps):\n",
" return 1./(np.exp(beta*eps) + 1.)\n",
"\n",
"def fermi_distribution_derivative(beta, eps):\n",
" return -beta/4. / np.cosh(0.5*beta*eps)**2\n",
"\n",
"def chi0_q0_integral(t, beta):\n",
"\n",
" def integrand(eps):\n",
" rho = get_density_of_states(eps, t)\n",
" df = fermi_distribution_derivative(beta, eps)\n",
" return -rho * df\n",
"\n",
" chi0_q0, err = quad(integrand, -2.*t, 2.*t)\n",
"\n",
" return chi0_q0 \n",
"\n",
"def find_Uc(t, beta):\n",
"\n",
" def root_function(U):\n",
" chi0_q0 = chi0_q0_integral(t, beta)\n",
" return 1 - U * chi0_q0\n",
"\n",
" Uc = brentq(root_function, 0, 100.)\n",
" return Uc"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = 1.0\n",
"beta = 5.0\n",
"U_vec = np.array([0., 1., 2., 3., 4., 5., 6., 7., 8.])\n",
"\n",
"for U in U_vec:\n",
" chi0_q0 = chi0_q0_integral(t, beta)\n",
" plt.plot(U, 1 - chi0_q0 * U, 'x')\n",
"\n",
"Uc = find_Uc(t, beta)\n",
"\n",
"plt.plot(Uc, 0., '+r', label=r'$U_c=%2.2f$' % Uc)\n",
"plt.ylabel(r'$1 - U \\chi_0(q=0)$')\n",
"plt.xlabel(r'$U$')\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = 1.0\n",
"beta = 5.0\n",
"n = 1.0\n",
"\n",
"U_vec = np.array([0., 1., 2., 3., 4., 5., 6., 7., 8.])\n",
"m_vec = np.linspace(-0.5, 0.5, num=100)\n",
"\n",
"for U in U_vec:\n",
" mu = 0.5*U \n",
" E_vec = np.array(\n",
" [ get_total_energy_mf_ref(t, beta, U, mu, n, m) \n",
" for m in m_vec ])\n",
" plt.plot(m_vec, E_vec, '-', lw=0.5, label=r'$U=%2.2f$' % U)\n",
"\n",
"mu = 0.5*Uc\n",
"E_vec = np.array(\n",
" [ get_total_energy_mf_ref(t, beta, Uc, mu, n, m) \n",
" for m in m_vec ])\n",
"plt.plot(m_vec, E_vec, 'r-', lw=1.5, label=r'$U_c=%2.2f$' % Uc)\n",
" \n",
"plt.xlabel(r'$m$')\n",
"plt.ylabel(r'$E[m]$')\n",
"plt.legend(loc='upper right')\n",
"plt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## TPRF tensor valued calculation"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"t = 1.0\n",
"#U = 6.1\n",
"U = 1.8\n",
"mu = 0.5*U\n",
"\n",
"m = 0.0\n",
"n = 1.0\n",
" \n",
"h_loc = -U*m*np.diag([+1., -1.]) + (0.5*U*n - mu) * np.eye(2)\n",
"T = - t * np.eye(2)\n",
"\n",
"from triqs_tprf.tight_binding import TBLattice\n",
"\n",
"t_r = TBLattice(\n",
" units = [(1, 0, 0)],\n",
" hopping = {\n",
" # nearest neighbour hopping -t\n",
" (0,): h_loc,\n",
" (+1,): T,\n",
" (-1,): T,\n",
" },\n",
" orbital_positions = [(0,0,0)] * 2,\n",
" orbital_names = ['up', 'do'],\n",
" )\n",
"\n",
"e_k = t_r.fourier(t_r.get_kmesh(n_k=(256, 1, 1)))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"k_vec, e_vec = np.array([ [k.value[0], e_k[k][0,0]] for k in e_k.mesh ]).T\n",
"plt.plot(k_vec, e_vec, '-')\n",
"plt.xlabel(r'$k$')\n",
"plt.ylabel(r'$\\epsilon(k)$');"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"╔╦╗╦═╗╦╔═╗ ╔═╗ ┌┬┐┌─┐┬─┐┌─┐\n",
" ║ ╠╦╝║║═╬╗╚═╗ │ ├─┘├┬┘├┤ \n",
" ╩ ╩╚═╩╚═╝╚╚═╝ ┴ ┴ ┴└─└ \n",
"Two-Particle Response Function tool-box \n",
"\n",
"beta = 5.0\n",
"nk = 256\n",
"nw = 800\n",
"norb = 2\n",
"\n",
"Approx. Memory Utilization: 0.06 GB\n",
"\n",
"--> fourier_wk_to_wr\n",
"--> fourier_wr_to_tr\n",
"--> chi0_w0r_from_grt_PH (bubble in tau & r)\n",
"--> chi_wk_from_chi_wr (r->k)\n",
"\n",
"chi0_q0 =\n",
"[[0.16215903 0. 0. 0. ]\n",
" [0. 0. 0.16215903 0. ]\n",
" [0. 0.16215903 0. 0. ]\n",
" [0. 0. 0. 0.16215903]]\n",
"\n",
"(0, 0, 0, 0) 0.16215902618552902\n",
"(0, 0, 0, 1) 0.0\n",
"(0, 0, 1, 0) 0.0\n",
"(0, 0, 1, 1) 0.0\n",
"(0, 1, 0, 0) 0.0\n",
"(0, 1, 0, 1) 0.0\n",
"(0, 1, 1, 0) 0.16215902618552902\n",
"(0, 1, 1, 1) 0.0\n",
"(1, 0, 0, 0) 0.0\n",
"(1, 0, 0, 1) 0.16215902618552902\n",
"(1, 0, 1, 0) 0.0\n",
"(1, 0, 1, 1) 0.0\n",
"(1, 1, 0, 0) 0.0\n",
"(1, 1, 0, 1) 0.0\n",
"(1, 1, 1, 0) 0.0\n",
"(1, 1, 1, 1) 0.16215902618552902\n",
"\n",
"chi0_q0 = 0.16215902618552902\n",
"chi0_q0_ref = 0.162159026185\n"
]
}
],
"source": [
"nw = 400\n",
"beta = 5.0\n",
"\n",
"from triqs.gf import MeshImFreq, Idx\n",
"wmesh = MeshImFreq(beta=beta, S='Fermion', n_max=nw)\n",
"\n",
"from triqs_tprf.lattice import lattice_dyson_g0_wk\n",
"g0_wk = lattice_dyson_g0_wk(mu=0., e_k=e_k, mesh=wmesh)\n",
"\n",
"from triqs_tprf.lattice_utils import imtime_bubble_chi0_wk\n",
"chi00_wk = imtime_bubble_chi0_wk(g0_wk, nw=1)\n",
"print()\n",
"print('chi0_q0 =\\n', chi00_wk[Idx(0), Idx(0, 0, 0)].real.reshape((4,4)))\n",
"\n",
"print()\n",
"import itertools\n",
"for idxs in itertools.product(list(range(2)), repeat=4):\n",
" print(idxs, chi00_wk[Idx(0), Idx(0, 0, 0)].real[idxs])\n",
" \n",
"chi0_q0_ref = chi0_q0_integral(t, beta)\n",
"\n",
"print()\n",
"print('chi0_q0 =', chi00_wk[Idx(0), Idx(0, 0, 0)][0,0,0,0].real)\n",
"print('chi0_q0_ref =', chi0_q0_ref)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"k_vec, chi0_vec = np.array([ [k.value[0], chi00_wk[Idx(0), k][0,0,0,0]] for k in e_k.mesh ]).T\n",
"plt.plot(k_vec, chi0_vec, '-')\n",
"plt.xlabel(r'$k$')\n",
"plt.ylabel(r'$\\chi_0(w=0, k)$');"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"H_int = 1.8*c_dag(0,0)*c_dag(0,1)*c(0,1)*c(0,0)\n",
"[1*c(0,0), 1*c(0,1)]\n",
"\n",
"U_abcd =\n",
"[[ 0. +0.j 0. +0.j 0. +0.j -1.8+0.j]\n",
" [ 0. +0.j 0. +0.j 1.8+0.j 0. +0.j]\n",
" [ 0. +0.j 1.8+0.j 0. +0.j 0. +0.j]\n",
" [-1.8+0.j 0. +0.j 0. +0.j 0. +0.j]]\n",
"\n",
"(0, 0, 0, 0) 0j\n",
"(0, 0, 0, 1) 0j\n",
"(0, 0, 1, 0) 0j\n",
"(0, 0, 1, 1) (-1.8+0j)\n",
"(0, 1, 0, 0) 0j\n",
"(0, 1, 0, 1) 0j\n",
"(0, 1, 1, 0) (1.8+0j)\n",
"(0, 1, 1, 1) 0j\n",
"(1, 0, 0, 0) 0j\n",
"(1, 0, 0, 1) (1.8+0j)\n",
"(1, 0, 1, 0) 0j\n",
"(1, 0, 1, 1) 0j\n",
"(1, 1, 0, 0) (-1.8+0j)\n",
"(1, 1, 0, 1) 0j\n",
"(1, 1, 1, 0) 0j\n",
"(1, 1, 1, 1) 0j\n"
]
}
],
"source": [
"gf_struct = [[0, [0, 1]]]\n",
"\n",
"from triqs.operators import n, c, c_dag, Operator, dagger\n",
"H_int = U * n(0, 0) * n(0, 1)\n",
"print('H_int =', H_int)\n",
"\n",
"from triqs_tprf.rpa_tensor import get_rpa_tensor\n",
"from triqs_tprf.rpa_tensor import fundamental_operators_from_gf_struct\n",
"\n",
"fundamental_operators = fundamental_operators_from_gf_struct(gf_struct)\n",
"print(fundamental_operators)\n",
"U_abcd = get_rpa_tensor(H_int, fundamental_operators)\n",
"\n",
"print()\n",
"print('U_abcd =\\n', U_abcd.reshape((4, 4)))\n",
"\n",
"print()\n",
"import itertools\n",
"for idxs in itertools.product(list(range(2)), repeat=4):\n",
" print(idxs, U_abcd[idxs])"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"chi_q0 =\n",
"[[ 0.17726126 0. 0. -0.05174012]\n",
" [ 0. 0. 0.22900138 0. ]\n",
" [ 0. 0.22900138 0. 0. ]\n",
" [-0.05174012 0. 0. 0.17726126]]\n",
"\n",
"(0, 0, 0, 0) 0.1772612564906984\n",
"(0, 0, 0, 1) 0.0\n",
"(0, 0, 1, 0) 0.0\n",
"(0, 0, 1, 1) -0.05174012291931889\n",
"(0, 1, 0, 0) 0.0\n",
"(0, 1, 0, 1) 0.0\n",
"(0, 1, 1, 0) 0.22900137941001733\n",
"(0, 1, 1, 1) 0.0\n",
"(1, 0, 0, 0) 0.0\n",
"(1, 0, 0, 1) 0.22900137941001733\n",
"(1, 0, 1, 0) 0.0\n",
"(1, 0, 1, 1) 0.0\n",
"(1, 1, 0, 0) -0.05174012291931889\n",
"(1, 1, 0, 1) 0.0\n",
"(1, 1, 1, 0) 0.0\n",
"(1, 1, 1, 1) 0.1772612564906984\n",
"\n",
"chi_SzSz_q0 = 0.11450068970500865\n",
"chi_q0_ref = 0.114500689705\n",
"\n",
"C/(1-UC) = 0.229001379409\n",
"C/(1-U^2C^2) = 0.17726125649\n",
"UC*C/(1-U^2C^2) = 0.0517401229191\n"
]
}
],
"source": [
"from triqs_tprf.lattice import solve_rpa_PH\n",
"chi_wk = solve_rpa_PH(chi00_wk, U_abcd)\n",
"\n",
"print('chi_q0 =\\n', chi_wk[Idx(0), Idx(0, 0, 0)].real.reshape((4,4)))\n",
"\n",
"print()\n",
"for idxs in itertools.product(list(range(2)), repeat=4):\n",
" print(idxs, chi_wk[Idx(0), Idx(0, 0, 0)][idxs].real)\n",
"\n",
"Sz = 0.5 * np.diag([+1., -1.])\n",
"chi_SzSz_wk = chi_wk[0,0,0,0].copy()\n",
"chi_SzSz_wk.data[:] = np.einsum('wkabcd,ab,cd->wk', chi_wk.data, Sz, Sz)\n",
"\n",
"print()\n",
"print('chi_SzSz_q0 =', chi_SzSz_wk[Idx(0), Idx(0, 0, 0)].real)\n",
"\n",
"# Eq. 7.43 Fazekas (additional 0.5 factor)\n",
"chi_q0_ref = 0.5 * chi0_q0_ref / (1. - U*chi0_q0_ref)\n",
"print('chi_q0_ref =', chi_q0_ref)\n",
"\n",
"C = chi0_q0_ref\n",
"\n",
"print()\n",
"print('C/(1-UC) =', C/(1.-U*C))\n",
"print('C/(1-U^2C^2) =', C/(1.-U*U*C*C))\n",
"print('UC*C/(1-U^2C^2) =', U*C*C/(1.-U*U*C*C))"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for idxs in itertools.product(list(range(2)), repeat=4):\n",
" k_vec, chi_vec = np.array([ [k.value[0], chi_wk[Idx(0), k][idxs]] for k in e_k.mesh ]).T\n",
" if np.max(np.abs(chi_vec)) == 0.:\n",
" continue\n",
" plt.plot(k_vec, chi_vec, '-', label=idxs.__repr__())\n",
"plt.legend()\n",
"plt.xlabel(r'$k$')\n",
"plt.ylabel(r'$\\chi(w=0, k)$');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{1}{1 + x} = \\sum_k (-x)^k\n",
"$$\n",
"\n",
"$$\n",
"e^x = \\sum_k \\frac{(-1)^k}{k!} x^k\n",
"$$\n",
"\n",
"$$\n",
"f(\\epsilon) = \\frac{1}{e^{\\beta \\epsilon} + 1}\n",
"= \\sum_k (-1)^k e^{k \\beta \\epsilon}\n",
"= \\sum_k (-1)^k \\sum_l \\frac{(-1)^l}{l!} (k \\beta \\epsilon)^l\n",
"$$\n",
"\n",
"$$\n",
"f( e \\mathbf{1} + \\eta \\sigma_x )\n",
"$$\n",
"\n",
"special case $[\\mathbf{1}, \\sigma_x] = 0$, anything commutes with $\\mathbf{1}$!\n",
"\n",
"Q: $[A, B] = 0$?\n",
"\n",
"Want\n",
"$$\n",
"\\lim_{\\eta \\rightarrow 0}\n",
"\\frac{ f( A + \\eta B ) - f(A) } { \\eta } ?=? B f'(A)\n",
"$$\n",
"in the general case $[A , B] \\ne 0$\n",
"\n",
"No what we want is slightly less general $A = A^\\dagger$ and $B = B^\\dagger$ and $B$ has only two non-zero terms that are unity.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.162159026185\n",
"0.162159026186\n"
]
}
],
"source": [
"def chi0_q0_integral_kspace(t, beta):\n",
" \n",
" def disp(k, t):\n",
" return -2.* t * np.cos(k)\n",
" \n",
" def integrand(k, t, beta):\n",
" ek = disp(k, t)\n",
" val = - fermi_distribution_derivative(beta, ek)\n",
" return val\n",
" \n",
" chi0_q0, err = quad(integrand, -np.pi, np.pi, args=(t, beta))\n",
" \n",
" chi0_q0 *= 1./(2.*np.pi)\n",
" \n",
" return chi0_q0\n",
"\n",
"\n",
"def fdm(E, beta):\n",
" from scipy.linalg import expm\n",
" I = np.eye(E.shape[0])\n",
" return np.mat(np.linalg.inv(expm(beta * E) + I))\n",
"\n",
"def dfdm(E, beta):\n",
" from scipy.linalg import coshm\n",
" B = np.mat(coshm(0.5*beta*E))\n",
" return -0.25 * beta * np.mat(np.linalg.inv(B*B))\n",
"\n",
"beta = 5.0\n",
"t = 1.0\n",
"\n",
"chi0_q0_ref = chi0_q0_integral(t, beta)\n",
"chi0_q0_ref2 = chi0_q0_integral_kspace(t, beta)\n",
"\n",
"print(chi0_q0_ref)\n",
"print(chi0_q0_ref2)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[]"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"